package com.swapper.math.ft;

import com.swapper.math.utils.BigDecimalUtils;

import java.util.Arrays;

@SuppressWarnings("all")
public final class FourierTransform {
  private static final int FLAG_DIRECT = 1;
  private static final int FLAG_INVERSE = -1;

  private static int checkLength(double[] realParts, double[] imaginaryParts) {
    int length = realParts.length;
    if (length != imaginaryParts.length) {
      throw new IllegalArgumentException("The length of the arrays(real and imaginary parts) is inconsistent.");
    }
    if (length < 1) {
      throw new IllegalArgumentException("The arrays(real and imaginary parts) is empty.");
    }
    return length;
  }

  public static String toArrayString(double[] realParts, double[] imaginaryParts) {
    int length = checkLength(realParts, imaginaryParts);
    String[] array = new String[length];
    for (int i = 0; i < length; ++i) {
      String crs = BigDecimalUtils.toPlainString(realParts[i]);
      String cis = BigDecimalUtils.toPlainString(imaginaryParts[i]);
      if (cis.startsWith("-")) {
        array[i] = crs + cis + "i";
      } else {
        array[i] = crs + "+" + crs + "i";
      }
    }
    return Arrays.deepToString(array);
  }

  public static void dft(double[] realParts, double[] imaginaryParts) {
    int length = checkLength(realParts, imaginaryParts);
    if (length > 1) {
      DFT(FLAG_DIRECT, length, realParts, imaginaryParts);
    }
  }

  public static void idft(double[] realParts, double[] imaginaryParts) {
    int length = checkLength(realParts, imaginaryParts);
    if (length > 1) {
      DFT(FLAG_INVERSE, length, realParts, imaginaryParts);
    }
  }

  public static void fft(double[] realParts, double[] imaginaryParts) {
    int length = checkLength(realParts, imaginaryParts);
    if (length > 1) {
      if ((length & (length - 1)) == 0) {
        if ((length & 0x55555555) != 0) {
          R4FFT(FLAG_DIRECT, length, realParts, imaginaryParts);
        } else {
          R2FFT(FLAG_DIRECT, length, realParts, imaginaryParts);
        }
      } else {
        DFT(FLAG_DIRECT, length, realParts, imaginaryParts);
      }
    }
  }

  public static void ifft(double[] realParts, double[] imaginaryParts) {
    int length = checkLength(realParts, imaginaryParts);
    if (length > 1) {
      if ((length & (length - 1)) == 0) {
        if ((length & 0x55555555) != 0) {
          R4FFT(FLAG_INVERSE, length, realParts, imaginaryParts);
        } else {
          R2FFT(FLAG_INVERSE, length, realParts, imaginaryParts);
        }
      } else {
        DFT(FLAG_INVERSE, length, realParts, imaginaryParts);
      }
    }
  }

  /**
   * The dft implement.
   *
   * @param F flag, {@link #FLAG_DIRECT} is fft, {@link #FLAG_INVERSE} is ifft.
   * @param L length.
   * @param X real part array.
   * @param Y imaginary part array.
   */
  private static void DFT(int F, int L, double[] X, double[] Y) {
    int I, J;
    double RAD, MR, MI;
    double[] X0 = X.clone();
    double[] Y0 = Y.clone();
    for (I = 0; I < L; ++I) {
      X[I] = 0D;
      Y[I] = 0D;
      for (J = 0; J < L; ++J) {
        RAD = -2 * Math.PI * F * I * J / L;
        MR = Math.cos(RAD);
        MI = Math.sin(RAD);
        X[I] += MR * X0[J] - MI * Y0[J];
        Y[I] += MR * Y0[J] + MI * X0[J];
      }
    }
    if (F == FLAG_INVERSE) {
      for (int i = 0; i < L; ++i) {
        X[i] /= L;
        Y[i] /= L;
      }
    }
  }

  /**
   * The radix-2 fft implement.
   *
   * @param F flag, {@link #FLAG_DIRECT} is fft, {@link #FLAG_INVERSE} is ifft.
   * @param L length, must be 2^n.
   * @param X real part array.
   * @param Y imaginary part array.
   */
  private static void R2FFT(int F, int L, double[] X, double[] Y) {
    double TEMP;
    int H, I, J, K, M, N;
    int[] REV = new int[L];
    for (I = 1, H = -1; I < L; I <<= 1) ++H;
    for (I = 0; I < L; ++I) {
      J = (REV[I >> 1] >> 1) | ((I & 1) << H);
      if (I < J) {
        TEMP = X[I];
        X[I] = X[J];
        X[J] = TEMP;
        TEMP = Y[I];
        Y[I] = Y[J];
        Y[J] = TEMP;
      }
      REV[I] = J;
    }
    double RAD, WNR, WNI, WR, WI, UR, UI, TR, TI;
    for (I = 2; I <= L; I <<= 1) {
      H = I >> 1;
      RAD = -2 * Math.PI * F / I;
      WNR = Math.cos(RAD);
      WNI = Math.sin(RAD);
      for (J = 0; J < L; J += I) {
        WR = 1D;
        WI = 0D;
        for (K = 0; K < H; ++K) {
          M = J + K;
          N = M + H;
          UR = X[M];
          UI = Y[M];
          TR = WR * X[N] - WI * Y[N];
          TI = WR * Y[N] + WI * X[N];
          X[M] = UR + TR;
          Y[M] = UI + TI;
          X[N] = UR - TR;
          Y[N] = UI - TI;
          TEMP = WR;
          WR = TEMP * WNR - WI * WNI;
          WI = TEMP * WNI + WI * WNR;
        }
      }
    }
    if (F == FLAG_INVERSE) {
      for (I = 0; I < L; ++I) {
        X[I] /= L;
        Y[I] /= L;
      }
    }
  }

  /**
   * The radix-4 fft implement.
   *
   * @param F flag, {@link #FLAG_DIRECT} is fft, {@link #FLAG_INVERSE} is ifft.
   * @param L length, must be 4^n.
   * @param X real part array.
   * @param Y imaginary part array.
   */
  private static void R4FFT(int F, int L, double[] X, double[] Y) {
    int J, K, M, N, Q;
    int I0, I1, I2, I3;
    double RAD, RAD1, RAD2, RAD3;
    double W1R, W2R, W3R, W1I, W2I, W3I;
    double F1R, F1I, F2R, F2I, F3R, F3I, F4R, F4I;
    for (M = 1, Q = L; M < L; M <<= 2) {
      N = Q;
      RAD = -2 * Math.PI * F / N;
      for (K = 0, Q >>= 2; K < Q; ++K) {
        RAD1 = K * RAD;
        RAD2 = RAD1 + RAD1;
        RAD3 = RAD1 + RAD2;
        W1R = Math.cos(RAD1);
        W1I = Math.sin(RAD1);
        W2R = Math.cos(RAD2);
        W2I = Math.sin(RAD2);
        W3R = Math.cos(RAD3);
        W3I = Math.sin(RAD3);
        for (I0 = K; I0 < L; I0 += N) {
          I1 = I0 + Q;
          I2 = I1 + Q;
          I3 = I2 + Q;
          F1R = X[I0] + X[I2];
          F3R = X[I0] - X[I2];
          F1I = Y[I0] + Y[I2];
          F3I = Y[I0] - Y[I2];
          F2R = X[I1] + X[I3];
          F4R = X[I1] - X[I3];
          F2I = Y[I1] + Y[I3];
          F4I = Y[I1] - Y[I3];
          X[I0] = F1R + F2R;
          F2R = F1R - F2R;
          F1R = F3R - F4I;
          F3R = F3R + F4I;
          Y[I0] = F1I + F2I;
          F2I = F1I - F2I;
          F1I = F3I + F4R;
          F3I = F3I - F4R;
          X[I2] = W2R * F2R - W2I * F2I;
          Y[I2] = W2R * F2I + W2I * F2R;
          if (F == FLAG_INVERSE) {
            X[I1] = W1R * F1R - W1I * F1I;
            Y[I1] = W1R * F1I + W1I * F1R;
            X[I3] = W3R * F3R - W3I * F3I;
            Y[I3] = W3R * F3I + W3I * F3R;
          } else {
            X[I1] = W1R * F3R - W1I * F3I;
            Y[I1] = W1R * F3I + W1I * F3R;
            X[I3] = W3R * F1R - W3I * F1I;
            Y[I3] = W3R * F1I + W3I * F1R;
          }
        }
      }
    }
    double TEMP;
    for (J = 1, K = 0, Q = L >> 2, N = L - 1; J < N; ++J) {
      M = Q;
      while (K >= 3 * M) {
        K -= 3 * M;
        M >>= 2;
      }
      K += M;
      if (J < K) {
        TEMP = X[J];
        X[J] = X[K];
        X[K] = TEMP;
        TEMP = Y[J];
        Y[J] = Y[K];
        Y[K] = TEMP;
      }
    }
    if (F == FLAG_INVERSE) {
      for (K = 0; K < L; ++K) {
        X[K] /= L;
        Y[K] /= L;
      }
    }
  }

  private static void SRFFT(int F, int L, double[] X, double[] Y) {

  }
}
